\documentclass[twoside,a4paper]{article}
\usepackage{geometry}
\geometry{margin=1.5cm, vmargin={0pt,1cm}}
\setlength{\topmargin}{-1cm}
\setlength{\paperheight}{29.7cm}
\setlength{\textheight}{25.3cm}

% useful packages.
\usepackage[UTF8]{ctex}	% 中文支持
\usepackage{listings} 	% 插入代码
\usepackage{xcolor} 	% 语法高亮
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{enumerate}
\usepackage{graphicx}
\usepackage{multicol}
\usepackage{fancyhdr}
\usepackage{layout}

% some common command
\newcommand{\dif}{\mathrm{d}}
\newcommand{\avg}[1]{\left\langle #1 \right\rangle}
\newcommand{\difFrac}[2]{\frac{\dif #1}{\dif #2}}
\newcommand{\pdfFrac}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\OFL}{\mathrm{OFL}}
\newcommand{\UFL}{\mathrm{UFL}}
\newcommand{\fl}{\mathrm{fl}}
\newcommand{\op}{\odot}
\newcommand{\Eabs}{E_{\mathrm{abs}}}
\newcommand{\Erel}{E_{\mathrm{rel}}}

\newcommand{\R}{\mathbb{R}}
\renewcommand{\proofname}{Proof.}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{exercise}{Exercise}[section]
\newtheorem{definition}{Definition}[section]

\begin{document}
	
	\pagestyle{fancy}
	\fancyhead{}
	\lhead{行一凡 (3190105815)}
	\chead{Project 3}
	\rhead{\today}
	
	\section{项目结构}
	程序所需头文件在 include 文件夹中，主程序文件在主目录下；三个文件夹 matlab、outputs、pictures 依次存放 matlab 绘图脚本、程序输出结果、生成图片。进入主目录下
	\begin{itemize}
		\item 输入 make story 编译生成 report.pdf 文件；
		\item 输入 make run 命令编译运行所有主程序，程序将结果输出到指定文件当中；
		\item 输入 make clean 清除可执行程序。
	\end{itemize}
	
	\section{程序结构}
	\begin{enumerate}[1.]
		\item include 文件中包含了程序所需头文件的相关代码。
		\begin{itemize}
			\item Vec.h 包装了 std::vector 为 Vec 类，用于储存向量，同时定义了所需的宏，以及变量别名。 TimeGrid 表示时间离散后的网格， GTimer 记录在网格中的位置，Func 和 Func2 定义了单变量向量函数和多变量向量函数;
			
			\item Factory.h 定义了生成单例类 TimeIntegrator 的对象工厂 GenericFactory ;
			
			\item Polynomial.h 定义了多项式类，用于求解多步法迭代的系数；
			
			\item EquationSolver.h 定义了求解线性方程组的类 LinearSolver 和执行牛顿迭代法以及获取 Jacobian 矩阵的类 FSolve ；
			
			\item TimeIntegrator.h 定义了单例类 TimeIntegrator 和其派生类 Multistep 和 RungeKutta ，它们又分别派生出具体的迭代方法；
			
			\item ODESolver.h 定义了具体的迭代过程，用户调用该类实现对求解方法的注册和调用方法求解 ODE ；
			
			\item Plotor.h 用于输出格式化的 matlab 代码，便于书写绘图程序。
		\end{itemize}
		
		\item 主目录中 Reportor.h 定义了测试方法的类，由各个测试程序直接调用来输出测试结果。
		\begin{itemize}
			\item Reportor.h 定义了用于测试的方法；
			
			\item test198.cpp 和 test199.cpp 分别输出两个初值条件下所有方法的测试结果到指定文件中；
			
			\item plot198.cpp 输出绘制 Euler 法，经典 RK 方法，Dormand-Prince 方法的绘图代码；
			
			\item winner.cpp 输出上述方法达到指定误差时的 CPU 时间的绘图代码。
		\end{itemize}
	\end{enumerate}
	
	\section{数值实现}
	对于多步法，程序使用多项式插值和积分方法来求解迭代系数，然后利用多步法的定义式定义隐函数
	\begin{equation}\label{key}
		F(u) = \alpha_su - k\beta_sf(u,t_n+k) + \sum_{i=0}^{s-1}\left[\alpha_iU^{n-i}-k\beta_if(U^{n-i},t_n-ik)\right]
	\end{equation}
	其中 $ (U^n,t_n) $ 是第 $ n $ 步迭代时的值和对应时间，使用中心差分方法计算 Jacobian 矩阵，然后用牛顿迭代法
	\begin{equation}\label{key}
		x_{n+1} = x_n - [J(x_n)]^{-1}F(x_n)
	\end{equation}
	求解 $ F(u)=0 $ 得到 $ U^{n+1} $ 。
	
	对于单步法，令 $ Y=(y_1,y_2,\cdots,y_s)^T $ ，同理利用定义式得到向量隐函数
	\begin{equation}\label{key}
		Y-F(Y) = 0
	\end{equation}
	同理计算 Jacobian 矩阵，然后应用牛顿迭代法求解得到 $ F(Y^*)=Y^* $ ，然后计算
	\begin{equation}\label{key}
		U^{n+1}=U^n+k\sum_{i=1}^{s}b_iy_i
	\end{equation}
	就得到下一步结果。
	
	\section{测试内容}
	由于系统中 $ u_3=u_6\equiv 0 $ ，因此在求解的 ODEs 系统中省略了对应的方程。
	
	\subsection{Test 1}
	使用初值
	\begin{equation}\label{key}
		U_0 = (0.994, 0, 0, -2.0015851063790825224)
	\end{equation} 
	时间步长 $ k=0.001 $ ，计算 $ T = 17.06521656015796 $ 的值，使用 $ U_0 $ 作为精确解，各个方法的结果如下
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/AdamsBashforth}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/AdamsBashforth198_T}
		\end{minipage}
		\caption{Adams-Bashforth}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/AdamsMoulton}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/AdamsMoulton198_T}
		\end{minipage}
		\caption{Adams-Moulton}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/BDFs}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/BDFs198_T}
		\end{minipage}
		\caption{BDFs}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/ClassicalRK}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/ClassicalRK198_T}
		\end{minipage}
		\caption{Classical RK}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/DormandPrince}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/DormandPrince198_T}
		\end{minipage}
		\caption{Dormand-Prince}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/ESDIRK}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/ESDIRK198_T}
		\end{minipage}
		\caption{ESDIRK}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/Fehlberg}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/Fehlberg198_T}
		\end{minipage}
		\caption{Fehlberg}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/GaussLegendre}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test198/GaussLegendre198_T}
		\end{minipage}
		\caption{Gauss-Legendre}
	\end{figure}
	
	% 换页
	\clearpage
	
	\subsection{Test 2}
	使用初值
	\begin{equation}\label{key}
		U_0 = (0.87978, 0, 0, -0.3797)
	\end{equation} 
	时间步长 $ k=0.004 $ ，计算 $ T = 19.14045706162071 $ 的值，使用外推法计算误差，各个方法的结果如下
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test199/ClassicalRK199_4}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test199/ClassicalRK199_T}
		\end{minipage}
		\caption{Classical RK}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test199/DormandPrince199_7}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test199/DormandPrince199_T}
		\end{minipage}
		\caption{Dormand-Prince}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test199/ESDIRK199_6}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test199/ESDIRK199_T}
		\end{minipage}
		\caption{ESDIRK}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test199/Fehlberg199_6}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test199/Fehlberg199_T}
		\end{minipage}
		\caption{Fehlberg}
	\end{figure}
	
	\begin{figure}[!htb]
		\centering
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test199/GaussLegendre199_3}
		\end{minipage}
		\begin{minipage}{0.45\linewidth}
			\centering
			\includegraphics[width=\linewidth]{pictures/test199/GaussLegendre199_T}
		\end{minipage}
		\caption{Gauss-Legendre}
	\end{figure}
	
	在 outputs 文件夹中记录了上面方法的平均收敛阶，简要整理如下
	\begin{center}
		\begin{tabular}{|c|c|c|c|c|}% 通过添加 | 来表示是否需要绘制竖线
			\hline  % 在表格最上方绘制横线
			方法 & Classical RK & Dormand-Prince & ESDIRK & Fehlberg \\
			\hline
			收敛阶 & 4.114355563 & 4.230412998 & 4.457673653 & 4.293852993\\
			\hline % 在表格最下方绘制横线
		\end{tabular}
	\end{center}
	
	\begin{center}
		\begin{tabular}{|c|c|c|c|}% 通过添加 | 来表示是否需要绘制竖线
			\hline  % 在表格最上方绘制横线
			方法 & Gauss-Legendre($ s=1 $) & Gauss-Legendre($ s=2 $) & Gauss-Legendre($ s=3 $) \\ 
			\hline
			收敛阶 & 2.001918705 & 3.96481235 & 1.787103494\\
			\hline % 在表格最下方绘制横线
		\end{tabular}
	\end{center}
	基本符合理论预计的收敛速度。其中 Gauss-Legendre 在 $ s=3 $ 时收敛很慢是因为其精度达到 $ 10^{-13} $ ，误差几乎不再变化。
	
	\clearpage
	\subsection{Test 3}
	测试 $ (10.198) $ 问题中 Euler 方法 24000 步，经典 RK 方法 6000 步，Dormand-Prince 方法使用自动控制步长 100 步，图像如下
	
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.5\linewidth]{pictures/plot198/Euler_plot198}
		\label{fig:eulerplot198}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.5\linewidth]{pictures/plot198/ClassicalRK_plot198}
		\label{fig:classicalrkplot198}
	\end{figure}
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.5\linewidth]{pictures/plot198/DormandPrince_plot100}
		\caption{100 步}
	\end{figure}
	
	% 换页
	\clearpage
	
	\subsection{Test 4}
	比较 max-norm 达到 $ 10^{-3} $ 时上述方法的 CPU 时间如下：
	
	\begin{figure}[!htb]
		\centering
		\includegraphics[width=0.7\linewidth]{pictures/winner/winner}
		\label{fig:winner}
	\end{figure}
	
	步长从 $ k=0.001 $ 开始减小，其中 Euler 方法虽然最快，但是一直不能达到精度要求；经典 RK 方法在 $ k=0.0005 $ 时达到精度要求，并且有较小的 CPU 时间；Dormand-Prince 方法误差平稳下降，但是 CPU 时间快速上升。因此经典 RK 方法在误差和 CPU 时间上都有更好的效果。
	
\end{document}